Traffic equilibrium cp PATH AMPL short i++3 */ param c {Nodes,Types}; /* delay coefficients for links in the network */ param gamma; /* gamma is a coupling coefficient. The delay on some of the links is increased gamma * delay on merging links */ var x {Nodes}; /* outside flow */ var z {Nodes}; /* inside flow */ var y {Nodes}; /* dual vars */ var fx {i in Nodes} = 1 + x[i] + x[i]^2; var fz {i in Nodes} = 1 + z[i] + z[i]^2; delay_o {i in Nodes}: 0 <= x[i] complements c[i,"oon"]*fx[i] + gamma * (1 + x[next(i,Nodes,3)] + x[next(i,Nodes,4)] + (x[next(i,Nodes,3)]+x[next(i,Nodes,4)])^2) + 10 * c[i,"ohy"]*(1 + x[i] + x[next(i,Nodes,3)] + x[next(i,Nodes,4)] + (x[i] + x[next(i,Nodes,3)] + x[next(i,Nodes,4)])^2) + 2 * gamma * fx[next(i,Nodes,3)] + c[next(i),"oby"] * (1 + x[i] + x[next(i,Nodes,4)] + (x[i] + x[next(i,Nodes,4)])^2) + 10 * c[next(i),"ohy"] * (1 + x[i] + x[next(i)] + x[next(i,Nodes,4)] + (x[i] + x[next(i)] + x[next(i,Nodes,4)])^2) + 2 * gamma * fx[next(i,Nodes,4)] + c[next(i,Nodes,2),"oby"] * (1 + x[i] + x[next(i)] + (x[i]+x[next(i)])^2) + 10 * c[next(i,Nodes,2),"ohy"] * (1 + x[i] + x[next(i)] + x[next(i,Nodes,2)] + (x[i] + x[next(i)] + x[next(i,Nodes,2)])^2) + 2 * gamma * fx[i] + c[next(i,Nodes,3),"ooff"] * fx[i] >= y[i]; delay_i {i in Nodes}: 0 <= z[i] complements c[i,"ion"]*fz[i] + gamma * fz[next(i)] + 10 * c[i,"ihy"]*(1 + z[i] + z[next(i)] + (z[i] + z[next(i)])^2) + 2 * gamma * fz[next(i)] + c[prev(i),"iby"]*fz[i] + 10 * c[prev(i),"ihy"] * (1 + z[i] + z[next(i,Nodes,4)] + (z[i]+z[next(i,Nodes,4)])^2) + 2 * gamma * fz[i] + c[prev(i,Nodes,2),"ioff"] * fz[i] >= y[i]; dem_cons {i in Nodes}: y[i] complements x[i] + z[i] = demand[i]; set InitPoints; param iptx {InitPoints}; param ipty {InitPoints}; param iptz {InitPoints}; ]]> set Types := ihy /* inside highway */ ion /* inside on-ramp */ ioff /* inside off-ramp */ iby /* inside bypass */ ohy /* outside highway */ oon /* outside on-ramp */ ooff /* outside off-ramp */ oby; /* outside bypass */ /* param demand := 1 1 2 8 3 1 4 8 5 1; demands are either ( .1, .2, .3, .4, .5) or ( 1, 8 1 8 1), with gamma set to 0, 0.5, or 4 */ param: Nodes: demand := 1 .1 2 .2 3 .3 4 .4 5 .5 ; param gamma := 0.5; param c default 1; param: InitPoints: iptx ipty iptz := run1 0 0 0 run2 0.5 0.5 0.5 run3 100 100 100 run4 0 0 0 run5 0 100 0 run6 0 1 0; for {run in InitPoints} { printf "Run %d:\n", run; let {i in Nodes} x[i] := iptx[run]; let {i in Nodes} y[i] := ipty[run]; let {i in Nodes} z[i] := iptz[run]; solve; display max{i in 1.._nccons} abs(_ccon[i]), min{i in 1.._ncons} _con[i].slack; }